from random import uniform

class ThrowMethod:
	def Throw(self):
		return 0.


class DirectPyThrow(ThrowMethod):
	def __init__(self, start = 0., end = 1.):
		self.start = start
		self.end = end

	def Throw(self):
		return uniform(self.start, self.end), uniform(self.start, self.end)

class DirectNaiveThrow(ThrowMethod):
	def __init__(self, seed = 43289, start = 0., end = 1.):
		self.start = start
		self.lens = end - start
		self.idum = seed
		self.m = 134456
		self.n = 8121
		self.k = 28411

	def Throw(self):
		self.idum = (self.idum*self.n + self.k) % self.m
		return self.start + self.lens*float(self.idum)/float(self.m), self.start + self.lens*float(self.idum)/float(self.m)


class MarkovThrow(ThrowMethod):
	def __init__(self, delta = .1, position = (1., 1.)):
		self.delta = delta
		self.x = position[0]
		self.y = position[1]

	def Throw(self):
		tmp_x = self.x + uniform(-self.delta, self.delta)
		tmp_y = self.y + uniform(-self.delta, self.delta)
		if (abs(tmp_x) < 1.) and (abs(tmp_y) < 1.):
			self.x = tmp_x
			self.y = tmp_y
		return self.x, self.y


##############################################################################
class Record:
	def __init__(self):
		self.hist = []

	def Clear(self):
		self.hist = []
	def Add(self, value):
		self.hist.append(value)

	def GetRecord(self):
		return self.hist


##############################################################################
class MonteCarlo_pi:
	def __init__(self):
		self.hits = 0
		self.rejects = 0
		self.x = 0.0
		self.y = 0.0

	def SetThrowMethod(self, throwMethod):
		self.throwMethod = throwMethod

	def Throw(self):
		self.x, self.y = self.throwMethod.Throw()

	def IsInCircleRange(self):
		return self.x**2. + self.y**2. < 1.

	def HitNotify(self):
		self.record.Add(1)

	def RejectNotify(self):
		self.record.Add(0)

	def SetRecord(self, record):
		self.record = record

	def GetRecord(self):
		return self.record.GetRecord()

	def Simulate(self, N):
		self.record.Clear()
		for i in range(N):
			self.Throw()
			if self.IsInCircleRange():
				self.HitNotify()
			else:
				self.RejectNotify()


class MonteCarlo_integral:
	def __init__(self):
		self.prob = 0.
		self.obser = 0.
		self.diff = 0.

	def SetFunction(self, prob, obser):
		self.prob = prob
		self.obser = obser
		self.diff = self.obser - self.prob

	def SetRecord(self, record):
		self.record = record

	def GetRecord(self):
		return self.record.GetRecord()

	def Throw(self):
		self.x = uniform(0, 1)**(1./(self.prob + 1.))

	def Simulate(self, N):
		self.record.Clear()
		for i in range(N):
			self.Throw()
			self.record.Add(self.x**self.diff)

##############################################################################
def Binning(base = 2, record = []):
#binning 2N array to N array
	tmp_sum = 0.0
	tmp_record = []
	count = 0
	for i in record:
		if count < base:
			tmp_sum += i
		else:
			tmp_record.append(tmp_sum/float(count))
			count = 0
			tmp_sum = i
		count += 1

	tmp_record.append(tmp_sum/float(count))
	return tmp_record	


def CalAvgStderr(record):
#calculate average stardard err
	sumAvg = 0.0
	sumSquar = 0.0
	count = 0
	for i in record:
		count += 1
		sumAvg += i
		sumSquar += i**2.
	avgStderr = ((sumSquar/float(count)-(sumAvg/float(count))**2.)/float(count))**(1./2.)
	return avgStderr


def CalStderrWithBinning(base, record):
	lens = 0
	avgStderrBinning = []
	numsBinning = []
	lens = len(record)
	if (lens > base):
		numsBinning.append(lens)
		avgStderrBinning.append(CalAvgStderr(record))
		while(lens > base):
			record = Binning(base, record)
			avgStderrBinning.append(CalAvgStderr(record))
			lens = len(record)
			numsBinning.append(lens)
	return numsBinning, avgStderrBinning



if __name__== '__main__':

	baseNum = 2
	N = baseNum**15
	delta = .3
	startPosition = (1., 1.)

	p = MonteCarlo_pi()
	#throwMethod = DirectPyThrow(0., 1.)
	#throwMethod = DirectNaiveThrow(12345, 0., 1.)
	throwMethod =  MarkovThrow(delta, startPosition)
	p.SetThrowMethod(throwMethod)
	p.SetRecord(Record())
	p.Simulate(N)
	numBinning, errBinning = CalStderrWithBinning(baseNum, p.GetRecord())
	for i in range(len(numBinning)):
		print numBinning[i], errBinning[i]
	import Gnuplot
	g = Gnuplot.Gnuplot()
	g.plot(errBinning)
	raw_input()

############################################################################
	#N = 1000
	#prob = 2.
	#obser = 3.

	#c = MonteCarlo_integral()
	#c.SetRecord(Record())
	#c.SetFunction(prob, obser)
	#c.Simulate(N)
	#print sum(c.GetRecord())/float(N)
